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Abstract 

A novel semi-Lagrangian method is introduced to solve numerically the Euler equation for ideal 
incompressible flow in arbitrary space dimension. It exploits the time-analyticity of fluid particle 
trajectories and requires, in principle, only limited spatial smoothness of the initial data. Efficient 
generation of high-order time-Taylor coefficients is made possible by a recurrence relation that fol¬ 
lows from the Cauchy invariants formulation of the Euler equation (Zheligovsky & Frisch, J. Fluid 
Mech. 2014, 749, 404-430). Truncated time-Taylor series of very high order allow the use of time 
steps vastly exceeding the Courant-Friedrichs-Lewy limit, without compromising the accuracy of 
the solution. Tests performed on the two-dimensional Euler equation indicate that the Cauchy- 
Lagrangian method is more — and occasionally much more — efficient and less prone to instability 
than Eulerian Runge-Kutta methods, and less prone to rapid growth of rounding errors than the 
high-order Eulerian time-Taylor algorithm. We also develop tools of analysis adapted to the Cauchy- 
Lagrangian method, such as the monitoring of the radius of convergence of the time-Taylor series. 
Certain other fluid equations can be handled similarly. 


1. Introduction 

As is well known, fluid flow can be characterised in terms of the current positions of fluid particles 
(Eulerian coordinates), or in terms of their initial positions (Lagrangian coordinates). For ideal 
(inviscid) incompressible fluid both formulations were introduced in the 18th century [15, 28]. From 
both a theoretical and a numerical point of view, the Eulerian formulation seems to have a significant 
edge because it gives an explicit quadratic expression for the time-derivative of the velocity; the 
Lagrangian formulation has even been qualified “an agony”, because of its complexity [36]. 

However, the Eulerian time-stepping methods have one serious well-known drawback, the Courant- 
Friedrichs-Lewy (CFL) condition [12], which constrains the time step to be less than a dimensionless 
constant multiplied by the time needed to sweep across the spatial mesh at the maximum flow speed. 
As a consequence, the complexity of computations with N collocation points in each spatial direc¬ 
tion is roughly O (IV 1 * * 4 ) in three dimensions. Hence, progress in high-resolution numerical simulation 
is creepingly slow, in spite of Moore’s law. (As we will see in Section 5, this is not the only drawback 
of Eulerian schemes.) 

Foremost because of the CFL constraint, which may conflict, for example, with the desire to 
quickly produce well-resolved numerical weather forecast, there has been a strong incentive to develop 
semi-Lagrangian (SL) schemes, in which some form of Lagrangian integration - with time steps not 
constrained by CFL - alternates with a reversion to an Eulerian grid (see, e.g., [40] and references 
therein). The SL algorithms used so far, e.g., in geophysical fluid dynamics, engineering, mechanics 
and plasma physics, were designed for situations where satisfactory results can be obtained with 
rather low-order temporal schemes. They are thus not appropriate for numerical experimentation on 
delicate questions, such as the issue of blow-up in three-dimensional flow (see, e.g., [22]). 




We propose a new SL algorithm, which we call the Cauchy-Lagrangian algorithm. It relies 
on Cauchy’s Lagrangian formulation of the equations of ideal incompressible flow [9, 20] and on 
recent results about the time-analyticity of Lagrangian trajectories [21, 44], The Cauchy-Lagrangian 
algorithm is particularly well-suited for problems where high precision is a prerequisite, and it is 
actually superior to Eulerian schemes. 

Cauchy’s 1815 Lagrangian equations are formulated in terms of the Lagrangian map from La¬ 
grangian to Eulerian coordinates, a —> x = x(a , t ), defined as the solution of the ordinary differential 
equation for fluid particle trajectory, x = v(x,t), with the initial condition x(a, 0) = a. Here, the 
dot denotes a (Lagrangian) time derivative. The three-dimensional Cauchy invariants equations are 
(see Section 2 for a simplified derivation): 

3 

^ W L x k x V L x k = u> (init \ det(V L cc) = 1, (1) 

k= 1 

where = V L x denotes the initial vorticity and V L is the Lagrangian gradient, i.e., the 

space derivatives in a. Since the r.h.s. of the first equation in (1) does not depend on time, its l.h.s. is 
obviously a Lagrangian invariant, whose scalar components are called the “Cauchy invariants”. These 
invariants were much later interpreted as a consequence of Noether’s theorem applied to a continuous 
symmetry of the Euler equation, called the relabelling invariance (see [20], Section 5.2). 

It was shown in [21, 44] that the Cauchy invariants equations (1) imply analyticity of fluid 
trajectories for initial conditions that have certain rather weak regularity. The proof of analyticity is 
derived from an explicit recurrence relation for coefficients of the time-Taylor series for the Lagrangian 
displacement $, = x — a. It can also be used to construct a numerical Lagrangian method of very 
high order in time in both two dimensions (2D) and three dimensions (3D). Its time steps are only 
constrained by the radius of convergence of the Taylor series (typically, the inverse of the largest 
initial velocity gradient). For more details, see Section 2 and [21, 44], 

The paper is organised as indicated hereafter. In Section 2 we recall some of the known results 
about Cauchy’s Lagrangian formulation and its application to the time-analyticity of the Lagrangian 
map. In Section 3 we describe the Cauchy-Lagrangian (CL) algorithm in detail: we begin with an 
overview and show that CL may be considered as a semi-Lagrangian algorithm of arbitrary high 
order (Section 3.1), present the interpolation technique for reverting to Eulerian coordinates at the 
end of each time step (Section 3.4) and show how to optimise the choices of the time step and of 
the order of the Taylor expansion truncation to minimise computational complexity (Section 3.5). 
Section 4 is devoted to testing the CL algorithm in the 2D case against various numerical methods, 
and to CPLI benchmarks. Section 5 is a comparison of high-order time-Taylor expansions in Eulerian 
and Lagrangian coordinates: the Lagrangian method suffers much less from the rounding errors. In 
Section 6 we determine how quickly we have to decrease the time step in the CL method because the 
radius of convergence R(t) of the time-Taylor series around time t generally shrinks as t increases, 
as more and more small-scale eddies (high spatial Fourier harmonics) are generated. Finally, in 
Section 7, we present concluding remarks and point out that the Cauchy-Lagrangian method is well 
adapted for certain other problems concerning ideal fluid flow. 

2. Mathematical background 

Here we recall some of the background material which is used to develop the Cauchy-Lagrangian 
method: derivation of the Cauchy invariants equations and the recurrence relation for time-Taylor 
coefficients from which follows the analyticity of the time-Taylor series (Section 2.1). We then 
present, in a heuristic way, a result allowing the determination of the radius of convergence of such 
Taylor series, an important new tool for analysing the output of Cauchy-Lagrangian computations 
(Section 2.2). 
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2.1. Cauchy’s Lagrangian formalism for ideal incompressible fluid flow 

The flow of an ideal incompressible fluid, described in Eulerian coordinates, is governed by the 
Euler equation 

d t v + (v ■ V)v = —Vp, V • v = 0, (2) 

where v(x, t ) is the velocity and p(x, t) the pressure (divided by the density, which in the incompress¬ 
ible case is constant). Here the spatial differentiation V is performed in the Eulerian coordinates x. 

The Eulerian equation (2) is notorious for difficulties in its investigation, both analytically and 
numerically. It turns out that its Lagrangian analogue can be once integrated in time, which yields 
the Cauchy invariants equations (1). Since these equations are central to our numerical method, we 
briefly recall how they are derived. 

In Lagrangian coordinates, fluid particles are characterised by their initial positions a and the time 
t. The subsequent positions of fluid particles are then described by the Lagrangian map a i —> x(a, t ). 
The spatial gradient in the Lagrangian coordinate a is denoted V L . 

Cauchy’s derivation of (1) is now briefly recalled [9]. Observe that the l.h.s. of (2) is the accel¬ 
eration x of the fluid particle; thus, x = —Vp. Multiplying this equation by the transpose of the 
Jacobian matrix V L x , we transform the Eulerian gradient in the r.h.s. into the Lagrangian one and 
obtain, in the component form: 

3 

^x^Xk = -Vj'p. (3) 

fc=i 

Cauchy then applies a Lagrangian curl to (3), and finds that the l.h.s. can be exactly integrated in 
time to yield the first equation of the set (1). The second one, the Jacobian equation, expresses the 
conservation of volume, equivalent to incompressibility. 

The Cauchy invariants equations are quadratically nonlinear equations, not resolved in the time 
derivatives; the Jacobian equation is also quadratic (in 2D) or cubic (in 3D). At first glance, it is 
unclear, in what respect (1) are advantageous compared to the Eulerian equations (2). To see what 
is gained, following [21, 44], we consider the equations for the displacement £ = x — a, that ensue 
from (1) and are in 3D as follows: 

3 

V L x £ + ^ V L 4 X V L e fc = ^ (init) , (4) 

fc =1 

Y L • £ • E (V^V^-V^V^) + det(V L O=0, (5) 

l<i<j<3 

and expand the displacement in the time-Taylor series 

OO 

= £ £<*>(«)«*. (6) 

s =1 


The structure of equations (4)-(5) enables us to derive recurrence relations for the coefficients ^ s \a), 
which in 3D take the following form [21, 44]: 

3 s —1 


V L x 4 s) 


v L • 4 s) 


-w y j v L zi m> x 

k=l m=1 


£ 

1<*<?<3 

£ 


s —1 

£ 


yLdm)yLt(s-m) 
v j Sj v i S j 


yLtWyL As-m) 
v i v j Sj 


m= 1 


( 7 ) 


( 8 ) 


i,j,k l-\-m-\-n=s 
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Here 8^ is the Kronecker delta and £ t] k is the totally antisymmetric tensor. In particular, for s — 1, 
(7)-(8) are equivalent to the relation = v- init k 

To determine $f s \ knowing its (Lagrangian) curl and divergence, is a Helmholtz-Hodge problem, 
with a unique solution for suitable boundary conditions. For example, for space-periodic flow we 
demand that the averages of all £^0 over the periodic box vanish, and construct the solution of 
(7)-(8) by taking the gradient of (8), subtracting from it the curl of (7) and thus obtaining a Poisson 
equation for £(*). This yields the following single (tensorial) recurrence relation for the (Lagrangian) 
gradients of the time-Taylor coefficients (s > 1) 


/ 


V^ s) = c, 

( 






E 

l</c<3, 
0 <m<s 


2 -^ (vy< m) ) 


( 9 ) 


■C, 


fiV 


\ 


Y. ((v J L d m) )v.yp m) - (vyytv^p’” 1 ) - Y 


l<i<i<3 
y 0<m<s 


l-\-m+n=s 


Here, C tl = V~ 2 V^V^ and denotes the second-order partial derivative d 2 /dctidctj. The oper¬ 

ator V~ 2 is the inverse of the (Lagrangian) Laplacian: given a periodic function / with zero spatial 
mean over the periodic cell, V~ 2 / is defined as the unique periodic function if with zero spatial 
mean, solving V 2 if = /. 

The explicit character of the recurrence relation (9) permits us to derive bounds for the timc- 
Taylor coefficients, provided the initial condition possesses a certain minimum smoothness, e.g., if 
the initial vorticity is Holder-continuous (i.e., belongs to the space C a for some 0 < a < 1) or 
belongs to the space of absolutely summable Fourier series (henceforth, ASFS). This way one can 
prove various analyticity theorems [21, 44], which state that the Lagrangian gradient of the Taylor 
series (6) converges in the space to which the initial vorticity belongs, provided |i| < inbound- Here, 
\t\ denotes the modulus of complex t and Abound a lower bound for the radius of analyticity , which 
is typically of the order of the inverse C a or ASFS norm of the initial vorticity. Within the disk of 
radius Abound, convergence of the time-Taylor series (6) and that of its spatial gradient is guaranteed 
and so is analyticity in time. The interested reader will find elementary proofs in [44], Of course, this 
is a key result for the construction of numerical solutions. Note that in 2D, by a theorem of Wolibner 
[42], when the initial vorticity is Holder-continuous, the vorticity stays smooth for all positive times. 
In 3D this problem is open. 

Suppose that, by using the Taylor series around t = 0 we have solved (1) in some real time interval 
[0,H], where t\ < Abound- If the solution at time t\ has the required smoothness, it can be used as 
a new initial condition for (1) at time t\. For this, Lagrangian coordinates must be reintroduced 
for fluid particles starting at t\. At t = H the new Lagrangian coordinates coincide, of course, 
with Eulerian coordinates. Hence we need to revert to Eulcrian coordinates at t,\. In principle, this 
can be done by composing the Lagrangian fields in the original coordinates with the inverse of the 
Lagrangian map, but there are other ways, as we will see. 

This procedure can be extended into a Weierstrass-type analytic continuation: one repeats the 
process at times t\ < as long as the successively constructed solutions have enough spatial 

smoothness to give time-Taylor series with a finite non-vanishing radius of convergence. In 2D, when 
the initial vorticity is Holder-continuous, the required smoothness will persist forever. In 3D loss of 
smoothness (blow-up) cannot be ruled out. 


2.2. The radius of convergence of the time-Taylor series 

When the initial vorticity is bounded in a suitable function space, such as C a or ASFS, we know 
that, for |i| < Abound, Ike gradient of the Taylor series (6) is guaranteed to converge in that space 
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for all a. Convergence may hold over a larger time interval, but without any proven control over 
the spatial smoothness of the Lagrangian map thus constructed. Knowing the radius of convergence 
of this Taylor series is important for numerical applications, even if |f| does not approach the edge 
of the disk of convergence: this parameter controls how fast the terms in the series tend to zero at 
large orders s. 

For a given starting point a, there exists a non-negative number R(a, 0), called the radius of 
convergence. It is the largest number such that the series (6) is absolutely convergent for any 
0 < \t | < R(a , 0). The second argument of the radius R indicates the time around which the Taylor 
expansion is performed. Note that the radius of convergence can be infinite. It is given by the 
Cauchy-Hadamard formula (see, e.g., [29]) 

* . =limsup|^(a)| 1/s = lim sup |£ (5) (a)| 1/s . (10) 

il(<X, 0) s—»oo S—K>o 5 > s 

It might be natural to define a radius of convergence i?(0) by taking the infimum of R(a , 0) over 
all initial points a. How can we measure this infimum in practice, given a set of Taylor coefficients 
^ s \a) up to a large value of s, and assuming that the a values are on a large set of collocation 
points? Following the usual numerical analyst’s approach, one might try to just evaluate the local 
radii of convergence R(a , 0) on a grid in the fluid domain and then find the minimum. This approach 
would have the advantage that it might indicate, what region in space is characterised by small values 
of the local radii of convergence, eventually perhaps leading to blow-up. However, the techniques 
for estimating the radius of convergence are numerically quite demanding, as discussed in Appendix 
B, suggesting that this approach is not numerically optimal. Furthermore, when the radius of 
convergence R(a , 0) is discontinuous, such an approach is clearly inappropriate. 

Even if £^(a) are analytic in a, we cannot guarantee the continuity of R(a, 0) (the series 
2~l a_a *l s t s , unrelated to the Euler equation, illustrates this difficulty). If R(a, 0) takes 
small values at exceptional off-grid points, we will miss them when measuring R( 0) numerically 
on a discrete grid. We find it therefore more suitable to employ the quantity 

^inf(O) = ess inf a R(a, 0), (11) 

where “ess” stands for “essential”, namely, small outlier values on a zero-measure set are disregarded. 

A simple way to evaluate i?i n f(0) makes use of the ordinary power series obtained by replacing 
£*A(a) by its L p norm, namely, 

OO 

£llc < ' ) (°)IU'- (i2) 

S=1 

We denote by R p the radius of convergence of the series (12). We assume that the initial vortic-ity 
is space-periodic and is either Holder-continuous or has an absolutely summable Fourier series. The 
same properties hold then for all the gradients of the Taylor coefficients V L £^(a). As a consequence, 
all the Taylor coefficients are in the Lebesgue spaces L p of functions, such that the pth power 

of their absolute value is integrable. In Appendix A we prove a theorem stating, roughly, that the 
radius -R in f(0) is also the radius of convergence of the series (12). More precisely, 

1. For any p > 1, we have i?i n f(0) > R p . 

2. If the displacement $,(a,t), defined by (6), is in L p for any 0 < t < -Rmf(O), then i?i n f(0) = R p . 
(The proof of this theorem does not involve the Euler equation.) 

For certain solutions of the Euler equation, the assumption in the second part of the theorem 
is known to be satisfied. For example, for 2D flow whose initial vorticity is Holder-continuous, 
Wolibner’s theorem [42] implies that the Lagrangian displacement is always in L p . In general, for 3D 
Euler flow, we do not know the spatial regularity properties of the displacement beyond the time of 
guaranteed analyticity and thus we only know that i? in f(0) > R p . 
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Figure 1: The successive disks of analyticity in the complex t plane. 

To gain heuristic insight into why the theorem holds, we observe that, since the radius of ana¬ 
lyticity depends only on the modulus of the Taylor coefficients, no generality is lost by taking them 
to be scalar and non-negative functions of a. For such a Taylor series, with radius of convergence 
R(a, 0), the coefficient of order s is roughly proportional to R(a,0)~ s (this would be exact for a 
geometric series). To calculate the L p norm of the coefficient of order s, we must integrate the pt\i 
power of this Taylor coefficient over a and take the pth root, ft seems plausible that, for large s, 
these integrals will be dominated by those a where R(a , 0) is close to its inffinum. We thus expect 
that the L p norm of the Taylor coefficient of order s is roughly proportional to (ess inf a A(a, 0)) _s , 
thereby yielding a radius of convergence given by R m {, namely, (11). 

3. Description of the Cauchy—Lagrangian algorithm 

In Section 3.1 we show that the basic ideas of the Cauchy-Lagrangian (CL) algorithm are quite 
simple, given our understanding of the analytical properties of the Lagrangian map, which were 
recalled in Section 2.1. Then, in Sections 3.2-3.5 we successively examine the constraints on the initial 
conditions (in connection with the spatial truncation errors), the temporal truncation errors, the 
reversion from the Lagrangian to the Eulerian representations and the minimisation of computational 
complexity. 

3.1. Overview: CL as a high-order semi-Lagrangian method 

Our purpose is to solve the Euler equation (2) with the initial condition u(:e, 0) = v( mit )(x) or 
uj(x,0) = cih mit )( 2 ;) in a time interval [0, T], In principle, we can just require that th im d or ub im b be 
in a function space that guarantees finite-time analyticity of the Lagrangian map, such as (for the 
vortic-ity) C a or the space of absolutely summable Fourier series ASFS. For practical reasons, more 
regularity is desirable (see Section 3.2). 

The simplest case is when T is less than the radius of analyticity bound Abound (0) of the tirne- 
Taylor series of the Lagrangian map, starting at t = 0. We then use the time-Taylor expansion of the 
displacement (6), together with the recurrence relation (9) to calculate as many Taylor coefficients 
as needed (see Sections 3.3 and 3.5). Here we just observe that, for spatial periodicity conditions, it 
is straightforward to calculate the Taylor coefficients using a pseudospectral method with dealiasing. 
In 3D, we can then obtain the vorticity in Lagrangian variables at time T by applying Cauchy’s 
vortic-ity formula [9]: 

u; L (a, T) = u/ init) (a) • V L a?(a, T). (13) 

(In 2D, this formula reduces to the constancy of the vorticity c o L (a, T ) = Cih im d(a).) Finally we can 
revert to Eulerian coordinates by composing this Lagrangian vorticity with the inverse Lagrangian 
map a{T) : x i-> a(x,T), which maps the fluid particle position x at time T to its initial position 
a , to obtain 

w(®,T) = w L (o(*,r),T). ( 14 ) 

Numerically, this is best done by interpolating from a Lagrangian to an Eulerian grid (see Section 3.4). 
Finally, the Eulerian velocity at time T , if needed, can be obtained either from the Eulerian vorticity 


6 



at time T (by inverting the curl, e.g., applying Fourier methods, or by interpolating to an Eulcrian grid 
the time derivative of the Lagrangian map, obtained by time differentiation of the Taylor expansion 
( 6 ))- 

Now, consider the case T > inbound(0) when convergence of the Taylor series is not guaranteed. 
Lack of convergence does not imply loss of regularity at some real time t between 0 and T. It may 
just be due to complex-time singularities whose moduli are less than or equal to T. Even if the 
Taylor series still converges, the sum might be insufficiently smooth to allow starting a new Taylor 
expansion (see end of Section 2.1). This is impossible in 2D with a Holder-continuous initial vorticity, 
but cannot be ruled out in 3D. The way to handle times T that are beyond the radius of analyticity 
bound inbound(0), is the multistep method, which essentially implements numerically the analytic 
continuation explained in Section 2.1. We pick a time t\ such that 0 < t\ < Abound (0) and start the 
process all over at time t\ using the Eulerian field v(x,ti) as the initial velocity (or u>(x,ti) as the 
initial vorticity). We then have, in the complex time plane, a new disk of guaranteed analyticity 
centred at t\ with radius inbound(U), which may well extend on the real positive time axis beyond 
the furthest point reached by the first disk of guaranteed analyticity. The process may be continued 
(see Fig. 1) with analyticity disks centred at various times 0 < U < t 2 < ..., as long as the solution 
remains in a space such as C“ or ASFS which guarantees that Abound > 0. 

The algorithm outlined above is illustrated by the flow chart in Fig. 2. 

To finish this overview, we observe that our Cauchy-Lagrangian 
method belongs to the family of semi-Lagrangian (SL) methods, of¬ 
ten employed, e.g., in meteorology, engineering, mechanics and plasma 
physics [5, 40]. An SL method solves an evolution equation along tra¬ 
jectories of particles for one time step and then reverts to the static 
regular Eulcrian grid, usually by interpolation. Unlike pure Eulcrian 
ones, SL methods are not subject to the CFL constraint [12]; conse¬ 
quently, their time steps are independent of the spatial resolution. For 
time-analytic Lagrangian fluid particle trajectories, the time step must 
just be smaller than the radius of convergence of the time-Taylor series, 
which is typically the inverse of the maximum velocity gradient. 

SL methods fall into two categories: upstream methods compute 
the so-called departure points, i.e., the positions of Lagrangian particles 
which are advected onto the regular grid during the time step, whereas 
for downstream SL methods the departure points of trajectories are 
on the regular grid, and the arrival points are on the distorted grid. 
The latter methods are simpler, because they do not require solving 
auxiliary problems to find departure points [30]. In this terminology, 
our CL method is downstream. Standard SL methods are typically of 
second-order accuracy in time; numerical schemes of higher orders are 
not widely used, being too complicated. Indeed, to advance by one 
time step, a high-order SL method of the traditional type would need 
information about the solution at many time levels; this would require 
many additional interpolations. In contrast, in the Cauchy-Lagrangian 
approach, the recurrence relation (9) allows us to construct very high- 
order schemes requiring only one interpolation per time step. 

3.2. Spatial truncation errors and the choice of initial conditions 
As we will see in Sections 4 and 5, our CL method is particularly well 
adapted and efficient for problems where high accuracy is important. 
Thus we want to keep all errors as (consistently) small as possible. 
This includes of course spatial truncation errors, stemming from the use of a finite number of modes. 


(^Initial Eulerian vorticityj 


Calculate Lagrangian 
time-Taylor coefficients 
at t tyI 


Estimate radius 
of convergence 
and choose time step 


Calculate Lagrangian 
map t n ->• t n 


Calculate Lagrangian 
vorticity at t = t n 


Transform to Eulerian 
vorticity at t — t n 



Figure 2: A flowchart of the 
CL algorithm operation. 


7 


















Actually, this is much easier if we limit ourselves to initial conditions that are analytic functions of 
the space variable. Then, the Eulerian solution remains analytic in space for at least a finite time 
(all times in 2D) [3]. During that time interval, the spatial Fourier transform of the solution falls 
off roughly as e _5 ^ fc , where k is the wave number and 5(t) is the radius of the tube of analyticity, 
the distance between the real spatial domain and the nearest complex-space singularity for the 
complexified spatial variables, ffence, truncation errors fall off then as e -5 W fcmax , where /c max is the 
maximum wave number present in the numerical truncation; thus they remain small as long as 5{t) 
is significantly larger than the spatial numerical mesh [41], 

In contrast, consider what happens when an initial condition has merely the minimum regularity 
required for ensuring the time-analyticity of the Lagrangian map, e.g., an initial vorticity is Holder- 
continuous. The spatial Fourier transforms of the velocity field during the time interval [0,T] will 
then fall off in a roughly algebraic way with the wave number k. Thus spatial truncation errors will 
fall off only algebraically with /c max . Achieving high accuracy may then require spatial resolutions 
that are beyond the reach of any existing computer. To avoid this, in our test computations we limit 
ourselves to spatially analytic initial data. 

One of the simplest space-analytic initial conditions that can be used for testing a CL code is 
the AB flow, discussed in Section 2.6.1 of [44], a steady-state solution of the 2D Euler equation. Its 
vorticity is = sin a; cosy. Although its Eulerian-coordinates dynamics is trivial, its Lagrangian 

dynamics is not completely trivial, but can still be determined analytically using elliptic integrals. 
More examples of flows with analytic initial data and truly non-trivial dynamics will be given in 
Section 4.1. 

3.3. Temporal truncation errors 

In numerical applications the Lagrangian displacement is approximated by the time-Taylor series 
(6) truncated to some finite order S: 


s 

£s(a, At) = £ (s \a)(At) s . (15) 

Which constraints on At and on S guarantee that the remainder of the series is sufficiently small? 

First, we demand that the series (6) be convergent, i.e., At should be less than the radius of 
convergence for all starting points a, and hence less than R\ n f, the inhmum over a of these radii 
(ignoring sets of zero measure). 

The determination of i?; n f proceeds as follows. Using the material of Section 2.2 and of Appendix 
A, we replace the full functional Taylor series by the series of the L p norms of the Taylor coefficients. 
This gives a series with radius of convergence R p . In 2D, R p = R inf and, in principle, in 3D we just 
have R p < R m f — however, since in 3D a blow-up of solutions of the Euler equation in the absence 
of solid boundaries perhaps never takes place, it is conceivable, at least for space-analytic initial 
conditions, that R p = R m {. The most convenient is to use the L 2 norm, which can be evaluated using 
ParsevaLs theorem. We thus have to determine the radius of convergence R 2 of a power series with 
positive coefficients. How to do this in practice is discussed in Appendix B. 

Second, we demand that the error due to truncation of the time-Taylor series be small. To see 
what this implies, we consider the remainder of the Taylor series and majorise its L 2 norm as follows: 

OO 

ll£(a»t) —£s(a,t)|| < Y, IIC (S, II(A*)', (16) 

s=S+l 

where || • || denotes the L 2 norm. The series in the r.h.s. has positive coefficients and its sum has 
a singularity at At = R 2 . Thus the sum is typically well approximated (at large orders and up to 
logarithmic factors) by a geometric series of ratio r = At/R 2 , and the L 2 norm of the remainder 



will typically be less than ||^' s )||(At)' s r/(l — r). For the values of the ratio r that minimise the 
complexity of the computation (see Section 3.5), r/( 1 — r) « 0.15. Hence, our criterion for keeping 
the truncation error small is 

||£ (S) ||(Af) s <e. (17) 

That is, we demand that, in the time-Taylor series, the last term kept before truncation be less, in L 2 
norm, than a prescribed accuracy e. Note that taking larger values of S allows us to advance forward 
with larger time steps. The choice of the optimal truncation order is discussed in Section 3.5. 


3.J h Reversion to an Eulerian grid by interpolation 


In our implementation of the Cauchy-Lagrangian method, the calculation of all the time-Taylor 
coefficients is done by pseudospectral techniques. For this to be done efficiently, the data at the 
beginning of a new time step must be known on a regular grid of N d collocation points, where N 
is the number of grid points per spatial period and d the dimension of space. This will be called 
here the uniform grid. Note that after the Lagrangian time-stepping is done, the new field, say, the 
vorticity, can either be determined in Lagrangian coordinates on the uniform grid or in Eulerian 
coordinates on a distorted grid, which is the image of the original grid by the Lagrangian map during 
the time interval At, shown in black in Fig. 3. Henceforth, we write just “image” for “image by the 
Lagrangian map”. 


To repeat the time-stepping process, we need to know the 
new Eulerian held, say, the vorticity, on a uniform Eulerian grid. 
We thus need to interpolate the Eulerian held from the distorted 
grid to the uniform grid. As we will see, the distortion of the grid 
remains moderate, so that the interpolation can be done using 
the cascade interpolation introduced in [37], which, in dimension 
d, is a superposition of d one-dimensional interpolations. We 
now briehy describe it. 

We begin with the 2D interpolation procedure, assuming 27r- 
periodicity and a time step from t — 0 to t — At. Positions of 
fluid particles at times 0 and At are a = (a, b ) and x = (x, y ) 
in the Lagrangian and Eulerian coordinates, respectively. The 
Lagrangian map takes a to x(a) = (x(a,b),y(a,b)). The uni¬ 
form grid is composed of the points a l3 = ( ai,bj ). Here, 
a* = 2ni/N and b 3 = 2nj/N, where the indices i and j run 
from 1 to the number of grid points per spatial period N. The 
same uniform grid is used for the Eulerian description after one time step. The N lines (a, b 3 ), where 
a varies continuously in the interval [0, 27r], are called the horizontal grid lines. Similarly, the N lines 
(a,;, b ) are called the vertical grid lines. 

The interpolation employs an intermediate hybrid grid, formed of intersections of the horizontal 
grid lines with the images of the vertical grid lines. This hybrid grid is shown by stars in Fig. 3. For 
the cascade interpolation, it is essential that each line of the former and latter type have a single 
intersection. This follows from a monotonicity result, discussed below. 

First, we determine the ^-coordinates of hybrid grid points. For each given pair of a horizontal 
and a vertical line this is done as follows. The images of the discrete grid points on the vertical 
line have coordinates (x n ,y n ). Along the image of the vertical line, in the neighbourhood of the 
intersection, a ID interpolation of x as a function of y is done; this gives the x value where y takes 
the value y 3 corresponding to the given horizontal line. Second, the vorticity, which is known at the 
same discrete locations (x n ,y n ) along the image of the given vertical line, is also ID-interpolated in 
y to the same values y 3 . Finally, along each horizontal line, the vorticity is ID-interpolated from 
the points on the hybrid grid to the points on the uniform grid. It is straightforward to extend the 
procedure to 3D with grid lines and grid planes and again only ID interpolations (see [37]). 



Figure 3: Distorted grid (black), image 
of the uniform grid (grey) by the La¬ 
grangian map, and hybrid grid (stars). 
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Now we turn to the issue of monotonicity of the map b i —> y(a j, b). First a qualitative observation: 
given the time-analyticity of the Lagrangian map and the fact that it is slightly smoother in space 
than C 1 , for small At the map is close to the identity map and its Jacobian matrix is close to 
the identity matrix. Hence, dy(a,b)/db > 0 for sufficiently small At, which implies the required 
monotonicity. This observation can be made quantitative by using the tools of [44] (Section 2.3). It 
is shown there (in 3D) that analyticity holds if TAt < T c = (8 — 5y / 2)/3, where T is the sum of the 
moduli of the Fourier coefficients of the initial vorticity. It is easy to show that the required positivity 
then holds at least to the same time. The positivity may cease to hold if At is just restricted to 
satisfy At < R m f, but in practice, we use values of At significantly smaller (see Section 3.5); in none 
of the simulations reported in Section 4 have we detected any lack of monotonicity. 

In the current version of the code, the ID interpolations are polynomial and use eight points. 
The interpolation procedure can be modified in a number of ways. We can interpolate first in 
the re-variable and then in the ^/-variable, and compare the results to estimate the accuracy of the 
interpolation. Instead of polynomial interpolation, we can use rational functions, splines, or any 
other ID interpolation algorithm. We can use different numbers of interpolation nodes in different 
regions of the periodicity square, depending, e.g., on how large is the vorticity gradient (note that in 
Fig. 7 of Section 4 the areas of large gradients are small for t > 3). For the interpolation we can use 
a mesh finer, than the one used to compute the Lagrangian map. Spatial derivatives of the vorticity 
can also be used for interpolation. We will experiment with these modifications in future studies. 

Errors of interpolations can be estimated on the fly as follows: the interpolation yields the 
vorticity on the target regular grid; from this, by FFT we can calculate its Fourier coefficients. By 
direct summation of the Fourier series (“slow Fourier transform”), we can then compute the vorticity 
on the distorted grid and find how closely it matches the one we employed for the interpolation. Such 
computations are CPU-intensive, but it suffices to perform them at those points where the largest 
interpolation errors are suspected, e.g., where high-order spatial derivatives of the vorticity are large. 

3.5. Choosing the optimal truncation order and time step 

In computations, the displacement is given by the time-Taylor series (15) truncated at order S. 
For a given accuracy £ , how should we choose S and the time step At to minimise the complexity 
of computations when integrating in a time interval [0,T]? We assume that the spatial resolution 
and thus the complexity of the spatial FFTs are fixed. The order of the interpolations needed for 
reversion to Eulerian coordinates is held fixed as well. Furthermore, we assume that the radius of 
convergence R- m f(t) of the time-Taylor series around time t does not vary significantly in the interval 
[0, T] and we denote it R. As we will see in Section 6, at least in two dimensions, R.mf(t) does not 
decrease in time dramatically, and hence this is a reasonable approximation. 

We prescribe the accuracy £ and apply the condition (17). The L 2 -norm series ||£^|| (At) s , 

which has positive coefficients and a radius of convergence R, can be approximated, at least for large 
orders, by the geometric series AY^^L 0 (At/R) s , where A is a positive constant. This exponential 
dependence on the order, 

||| (s) || « AR~ S , (18) 

is well supported by the numerical results on high-order time-Taylor expansions of Section 5 (see, in 
particular, Fig. 11). 

For large truncation orders S, given the nonlinearity of order d stemming from the determinant in 
(5), the number of operations is 0(5^). Here d is the space dimension. The number of time steps for 
integration to time T is T / At, (17) and (18) yield At = R(e / A) l ^ s . Hence, the number of operations 
needed to reach a given time is proportional to S d (e / A)~ l A. Minimising this expression over S, we 
hnd the optimal order S ~ —(1/d) In (e/A) and At = Re~ d . For small S and large N, the complexity 
SN d In N of FFTs needed for applying the Calderon-Zygmund operators involved in the expression 
of becomes dominant. By similar arguments, the optimal order is then S ~ — In (e/A). When 
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Figure 4: Temporal variation of the time step of the CL code for the 4-mode (a) and random initial conditions (b). 


both complexities are essential, the optimal order of expansion, S , satisfies the inequality 

— (1/d) hx(e/A) < S < — In (e/A). (19) 

4. Testing the Cauchy—Lagrangian numerical method in two dimensions: 

Comparison with Eulerian simulations 

In Section 4.1 we specify the test flows and describe three Eulerian algorithms that have been used 
for comparison with the Cauchy-Lagrangian algorithm. Validation, with emphasis on accuracy, is 
presented in Section 4.2. Efficiency of the CL algorithm is discussed in Section 4.3. In Section 4.4 we 
discuss spatial truncation artefacts. All computations reported in this section arc in double precision. 

4-1. Flows and numerical methods used for validation 

For the reason explained in Section 3.2, all our tests of the CL method have been done using 
2D flows with analytic initial data having non-trivial dynamics in both Eulerian and Lagrangian 
coordinates. Two different flows, called the test flows, were used as initial condition. The first one 
is a very simple flow, here called the “4-mode” flow, with the initial vorticity 

_ CQS x _|_ CQS y _|_ q g cos 2 X + 0.2 cos 3x. (20) 

The second one is a particular realisation of the so-called “random” flow [45], used in Section II. C 
of [39], where the time evolution of the latter flow is presented in detail. The characterisation of the 
random flow is best done in the Fourier space, where each wave vector consists of a pair of signed 
integers k = {k\, k- 2 ), conventionally combined into shells K < |fc| < K + 1 where K is integer. 
Each such shell involves N(K) Fourier harmonics. For all k in the K th shell, the Fourier coefficients 
cDfc of the initial vorticity are assigned the same modulus 2K‘I 2 exp(— K 2 /A)/N(K) and phases that 
are uniformly and independently distributed in the interval [0, 27 t], except that the coefficients for 
opposite wave vectors are complex conjugate, so that the flow is real. 

We have used three programs to solve the 2D Euler equation, realising the CL algorithm with 
the Taylor series truncated to order S (denoted CLA), described in Section 3, and three algorithms 
working exclusively in Eulerian variables, namely, the Runge-Kutta algorithms of orders two and 
four (denoted RK2 and RK4, respectively) and the algorithm relying on the Eulerian time-Taylor 
expansion truncated to order S (denoted ETS 1 ). The CL algorithms used in this section are CL8, 
CL16 and CL24 with time steps, chosen as explained in Section 3.5, and with the accuracy in (17) 
set to £ = ICC 12 . Figure 4 shows the evolution of the time steps for the test flows. In principle, they 
are allowed to vary in time, but, as seen, they change very moderately. This is apparently related to 
the very slow temporal variations of the radii of convergence of the timc-Taylor series, discussed in 
Section 6. 
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Our calculations have mostly been performed with spatial resolutions of 1024 2 and 2048 2 Fourier 
harmonics before dealiasing (performed using the standard 2/3-rule applying to quadratic problems). 
One calculation had a resolution of 8192 2 harmonics. 

Concerning the standard RK2 and RK4 algorithms, we only comment on the values of the time 
steps, taken constant in our calculations. For RK2, the smallest time step that we were able to 
run in reasonable CPU time was At = 10 -4 . As we will see, this has yielded accuracies in the 
range 10 -6 — 1CU 9 , significantly worse than the ones obtained with the other algorithms, requiring 
comparable or smaller CPU time. For RK4, the time steps used are listed in Table 1. They have 
been chosen as the largest value that gives stable integration, so that doubling of the time step causes 
numerical overflow. The corresponding Courant numbers (based on the maximum velocity and the 
maximum wave number after dealiasing) are indeed slightly in excess of unity. 

The Eulerian Taylor (ET) algorithm deserves more explanations. By analogy with time integra¬ 
tion exploiting analyticity of Lagrangian fluid particle trajectories, one can perform Eulerian time 
integration exploiting the analyticity in time of the velocity and the vorticity when the initial data 
are analytic in space and remain so for at least a finite time. The ET algorithm was actually applied 
long ago to study real and complex time singularities of the Taylor-Green 3D vortex flow [6]. This 
was done with a single time step from 0 to t, followed by attempted analytic continuation in t beyond 
the disk of convergence using, e.g., Pade approximants. Unfortunately, this method failed to give 
clear evidence for or against finite-time blow-up in 3D. Here, we are interested in the 2D multistep 
form of the ET algorithm applied to the vorticity formulation of the Euler equation 

Did 

— + (v • V) u; = 0, \/xv = coe z , (21) 

where u denotes the only component of the vorticity and e z is the unit vector in the z-direction. We 
assume that the vorticity held to is known at time t, and Taylor-expand tu(t + At) in powers of At: 

OO 

u(t + At) = u s (At) s . (22) 

s=0 

After substitution in (21), we obtain the following recurrence relations (where ojq = oo(t)): 

S 

(s + 1) CU S+ 1 + (v m • V) Vs-m = 0, Vxv m =u m e z , s = 0,1,.... (23) 

m =0 

These recurrence relations look significantly simpler than the Lagrangian recursion relation (9). 
Actually, the Lagrangian one is much better conditioned than the Eulerian ones, as we shall see in 
Section 5. 

Time-stepping is done, as in the CL method, by applying the Taylor series truncated to a finite 
order S and choosing an appropriate time step At. Here, we report results obtained with an 8-term 
truncation, called ET8. Higher-order truncations have special problems, that will be discussed in 
Section 5. 

Given the analyticity in time of the Eulerian solution (for all times in 2D and a finite time in 3D 
[3]), we certainly must require that the time step be smaller than the (Eulerian) radius of convergence 


Table 1: Time step At for RK4 and ET8 algorithms. 


Resolution 

1024 2 2048 2 

4-mode initial how 

random initial how 

0.0025 0.00125 

0.002 0.001 
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of the Taylor series. However, no Eulerian analogue of the condition (17) is available, and much 
smaller time steps must be chosen; typically Courant numbers of order unity are needed. The 
reason, however, is not stability, as in low-order schemes, but the rounding errors, as explained in 
Section 5. In practice, for ET, we use the time steps of the same length as for RK4, given in Table 1. 

Finally, a few words about spatial truncation errors. Since we restrict ourselves to analytic initial 
data, the solutions are also analytic in both Lagrangian and Eulerian coordinates. The Fourier 
transforms of such functions decrease exponentially with the wave number k as e~ 5k , where S is the 
radius of the cylinder of analyticity (also sometimes called “width of the analyticity strip”). As 
shown in [41], 6 can be estimated from the high- K asymptotics of the vorticity spectrum (also called 
“enstrophy spectrum”), 

E U (K) = \ Y. I a *l 2 - < 24 ) 

I<<\k\<K+l 

where cD*. are the spatial Fourier coefficients of the vorticity c o. The large-wave-number behaviour of 
the vorticity spectrum is typically 

E U (K) ~ K a e~ 2SK . (25) 

Here, the distance 5 tells us how far away from the real domain the nearest complex singularities 
are, while the exponent a contains information about the type of singularities [35]. As a consequence 
of (25), the relative spatial truncation errors are typically estimated as e -5fcmax , where k max is the 
maximum wave number after dealiasing, k max = At/3, where N is the number of grid points per spatial 
period. The numerical determination of 5 from the vorticity spectrum is done by the same fitting 
method as for the numerical determination of radii of convergence (see Appendix B). When using 
the Cauchy-Lagrangian method, at each new time step, the Lagrangian and Eulerian descriptions 
coincide and, as a consequence, the Lagrangian and Eulerian 5 are close. 

4-2. Validation of the CL algorithm and accuracy of agreement 

For each test flow (see Section 4.1), we have to determine up to what maximum time £g n we can 
integrate the 2D Euler equation without encountering excessive loss of accuracy. This will of course 
depend on the resolution and on how rapidly the Fourier coefficients fall off with the wave number k. 
As explained at the end of Section 4.1, this fall-off is controlled by the radius of spatial analyticity, 
measuring the distance S(t) of the nearest complex-space singularities to the real-space domain. To 
monitor S(t), we have used the Cauchy-Lagrangian method to integrate the Euler equation with 
the resolution of 8192 2 Fourier harmonics. The measured 5(t) for the two test flows are shown in 
Fig. 5. Each point on these graphs is obtained by processing the vorticity spectrum E UJ (K) at the 
corresponding time by the fitting technique of Appendix B. Examples of vorticity spectra are shown 
in Fig. 6. 

We determine from Fig. 5 that, at the largest times shown for both flows, 5 is about ICE 2 . Since 
fcmax = 8192/3 ~ 2731, the resulting relative truncation error is about e _<5 ^ fcmax « 2 x ICE 12 . (The 
product S(t)k max is sufficiently large to use this asymptotic wave number formula.) We also note that 
for the lowest resolution used in our simulations, namely, 1024 2 harmonics, the 4-mode flow achieves 
comparably small spatial truncation errors only up to t — 2. If, however, we just request a level of 
spatial truncation error of about ICE 4 , which is not visually detectable on contours of the vorticity 
and of its Laplacian, we can extend the computations to about t — 4. We will come back to this in 
Section 4.4, when discussing truncation artefacts. 

Before comparing the results for the different methods, it is of interest to display the evolution of 
the test flows graphically. For the random flow, this was done in Section II.C of [39]. Here we focus 
on the 4-mode flow. 

The evolutions of the vorticity and of its Laplacian for the 4-mode initial flow are shown in Figs. 7 
and 8, respectively. We also performed high-resolution (8192 2 ) simulations for the same initial flow. 
The flow organises itself into a periodic array of large eddies centred at the centre and corners of 
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Figure 5: Radius 6(t) of the analyticity cylinder for the 4-nrode (a) and the random (b) initial condition. CL8 
algorithm. Resolution: 8192 2 harmonics. 


(a) 
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Figure 6: Spectra of the vorticity for the 4-mode initial condition at several times, as labelled. CL8 algorithm. 
Resolution: 8192 2 harmonics. 

the displayed periodicity box. Fine structures involving large vorticity gradients develop at the 
boundaries between these large eddies. 

At t = 2, we see in plots of vorticity isolines, and more clearly in the isolines of its Laplacian, 
two thin vortical sheets (related by the central symmetry of the flow), which are roughly inclined at 
+45° with respect to the x-axis and two other ones, roughly inclined at -45°, which are not quite as 
fine. Around t = 3.7, the former and the latter have thinned to comparable scales. The thickness of 
such structures can be interpreted as a rough measure of the distance to complex-space singularities, 
because it determines the strength of spatial derivatives in the real-space domain. So, here we 
have a pair of singularities (actually, singular manifolds) at complex-conjugate spatial locations that 
competes with another similar pair, with eventually the one pair furthest from the real domain 
catching up with the one which was closer. This singularity catch-up explains the bending of the 
plot in Fig. 5 around t = 3.7. Another feature shown in Fig. 8(b), after the singularity catch-up, 
are oscillations superimposed on a roughly exponential energy spectrum. Detailed examination of 
contours of the Laplacian of vorticity at the highest resolution (8192 2 , not shown) reveals that the 
thin sheets seen in Fig. 8, which are inclined at -45° have a doublet structure, while those inclined 
at +45° have a singleton structure. The doublet structure gives rise to oscillatory interference in the 
Fourier transform. 

We now show that the CL results agree with those obtained by traditional methods for both 
test flows within the time intervals where spatial truncations effects are negligibly small. Table 2 
shows the vorticity discrepancies between results of different algorithms, namely, the maximum over 
space of the difference of the vorticity calculated with two different codes. This is done for the two 
test flows at different spatial resolutions and different output times. We also performed energy and 
enstrophy conservation tests, whose results are presented in Tables 3 and 4, respectively. (We will 
return to conservation issues in Section 4.4.) The main result is that the vorticity discrepancies 
between the CL results and those obtained by the standard Eulerian RK4 method are very small, 
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Figure 7: Isolines of the vorticity (step 0.5) computed by CL8 with the resolution of 1024 2 harmonics for the 4-mode 
initial condition. 



Figure 8: Isolines of the Laplacian of the vorticity (step 1 for t = 0; 2 for t = 1; 10 for t = 2; 60 for t = 3; 500 for 
t = 4 and 4.1), computed by CL8 with the resolution of 1024 2 harmonics for the 4-mode initial condition. 


around 10 -1 °-10~ 12 , except for later times, when spatial truncation errors become important. Note 
that the discrepancy with the RK4 results is roughly the same for CL8, CL16 and CL24. Choosing 
between these methods is thus only a matter of efficiency, as discussed in the next section. 

4-3. Efficiency of the CL method 

Before discussing the relative performances of the various algorithms, let us note that most of 
the computations were made on a computer - marketed in 2011 - that has an Intel Core i7-3960X 
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Table 2: Vorticity discrepancies between solutions (maximum over the periodicity box of the absolute value of the 
difference), computed by different methods for the test flows at various times and resolutions. 


Run 

t 

CL8-RK4 

CL16-RK4 

CL24-RK4 

ET8-RK4 

RK2-RK4 

4-mode, 

1 

1.5-10 -12 

1.5-10~ 12 

1.5-HT 12 

1.5-10~ 12 

10” 9 

1024 2 

3 

5-10 -11 

2-10~ n 

3-10” 11 

5-10- 12 

1.5-HT 9 


5 

6-10 -5 

5-10 -5 

5-10 -5 

2-10 _1 ° 

3-10" 9 

4-mode, 

1 

3-10" 12 

5-10” 12 

8-10 -12 

2-10~ 12 

10” 8 

2048 2 

3 

1.5-10 -10 

7-10" 11 

o 

1 

o 

3-10 -11 

o 

1 

00 


5 

8-10~ 4 

2-10 -4 

4-10" 4 

4-10" 8 

4-10" 6 

random, 

0.2 

3-10~ 10 

2-10 _1 ° 

2-10” 10 

2-10” 10 

3-10" 7 

2048 2 

0.6 

7-10 -8 

2-10 -8 

2-10” 9 

2-10" 9 

10~ 6 


1 

2-10 -3 

2-10 -3 

2-10 -3 

2-10~ 6 

4-10" 5 


Table 3: Energy conservation errors. 


Run 

t 

CL8 

CL16 

CL24 

RK2 

RK4 

ET8 

4-mode, 

1 

-3-10" 14 

io- 15 

< 10~ 15 

io- 11 

2-10 -14 

< 10“ 15 

1024 2 

3 

-2-10" 14 

7-10" 14 

2-10~ 13 

-2-10" 11 

4-10~ 14 

< IO" 15 


5 

-5-10" 9 

-IO” 8 

1 

o 

1 

00 

-4-10” 11 

7-10~ 14 

2-10~ 15 

4-mode, 

1 

-3-10" 14 

10- 15 

10~ 15 

10~ n 

10“ 15 

< 10~ 15 

2048 2 

3 

-2-10" 14 

10- 15 

10- 15 

—2-10 -11 

< 10~ 15 

-2-10' 15 


5 

—6-10 -11 

1 

o 

l 

o 

1 

o 

1 

o 

_ 4 .ll - 11 

10~ 15 

—5-10” 15 

random, 

0.2 

KBS 

-2-10" 15 

-2-10" 15 

—7-10” 11 

IO" 14 

-2-10' 15 

2048 2 

0.6 

-4-10" 14 

-4-10" 14 

—3-10” 14 

2-10” 11 

3-10~ 14 

-3-10" 14 


1 


IO" 10 

-2-10" 10 

6 -11 —10 

-3-10" 13 

—3-10~ 13 


Table 4: Enstrophy conservation errors. 


Run 

t 

CL8 

CL16 

CL24 

RK2 

RK4 

ET8 

4 - mode , 

1 

- 4 - 10" 14 

< 10~ 15 

< 10“ 15 

3 - 10” 11 

2 - 10 -14 

< 10~ 15 

1024 2 

3 

- 2 - 10" 12 

-io - 12 

— 3 - 10~ 12 

6 - 10” 11 

3 - 10~ 15 

— 3 - 10 -15 


5 

- 2 - 10" 6 

-IO" 6 

- 2 - 10” 6 

4 - 10” 10 

— 2 - 10 -12 

- 5 - 10' 15 

4 - mode , 

1 

- 4 - 10" 14 

< 10" 15 

< 10~ 15 

3 - 10” 11 

io - 15 

< 10“ 15 

2048 2 

3 

- 3 - 10" 14 

- 8 - 10" 14 

— 10 -14 

O 

o 

1 

5 - 10 -15 

- 6 - 10' 15 


5 

- 2 - 10" 8 

— 5 - 10 -9 

- 7 - 10' 9 

4 - 10” 10 

10~ 13 

- 4 - 10' 14 

random , 

0.2 

io - 14 

- 2 - 10" 14 

- 2 - 10” 14 

— 3 - 10~ 10 

10~ 13 

- 2 - 10” 14 

2048 2 

0.6 

— 6 - 10 -13 

- 2 - 10" 13 

- 2 - 10' 13 

4 - 10 -9 

10- 13 

- 2 - 10" 13 


1 

-IO” 8 

- 10“ 8 

2 - 10 -9 

3 - 10“ 9 

— 8 - 10~ 12 

- 2 - 10" 12 
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Table 5: CPU times needed to reach various output times for the codes applied to the test flows. 


Run 

t 

CL 8 

CL16 

CL24 

RK2 

RK4 

ET 8 

4-mode, 

1 

lulls 

27s 

32s 

92m50s 

7m33s 

16m22s 

1024 2 

3 

2m56s 

lml5s 

lml5s 

278m20s 

22m26s 

48ml s 


5 

4m55s 

2 m 12 s 

2m08s 

464m0s 

37m22s 

81m50s 

4-mode, 

1 

4m20s 

lm59s 

2ml7s 

402m20s 

65m01s 

140ml Os 

2048 2 

3 

12m30s 

5m24s 

5ml8s 

1206m50s 

194m20s 

420m20s 


5 

20m56s 

9m23s 

9m07s 

2011 m 20 s 

323m40s 

700ml0s 

random, 

0.2 

3m05s 

lm29s 

lm32s 

81m40s 

16ml9s 

34m53s 

2048 2 

0.6 

8m42s 

3m57s 

4m32s 

244m50s 

48m20s 

103m40s 


1 

13m48s 

6m25s 

6m51s 

408m20s 

81ml0s 

163m40s 


processor (64 bit data width, frequency 3.3 GHz up to 3.9 GHz in the Turbo regime, 15 Mbyte L3 
cache). It has 6 cores with hyper-threading, which does not actually matter for our application, since 
we have so far experimented only with sequential versions of the codes. Parallelising such a code is 
standard - at least for the Lagrangian time-stepping phase - and requires the same tools as for the 
parallelisation of any spectral code. 

The CPU timings for the integration to various output times for our codes are shown in Table 5. 
The timings for the Lagrangian CL codes are way below those for the best Eulerian code: at the 
resolution of 2048 2 harmonics, CL24 needs at least an order of magnitude less CPLI time than the 
commonly used RK4; the ET 8 algorithm is at least 2 times slower than RK4, and RK2 is much 
slower. 

We also observe that CL 16 and CL24 have comparable efficiency and are significantly faster than 
CL 8 . Hence the optimal number of terms in the time-Taylor series is around S = 20. This is 
consistent with the arguments of Section 3.5 that give bounds for the optimal order, involving an 
unknown constant A. 

All the results given so far are for the choice £ = 10 ~ 12 for the accuracy parameter that controls 
the truncation errors on the time-Taylor series (see (17)). It is also of interest to measure the CPLI 
time dependence on e. For the four-mode initial conditions at resolution 1024 2 integrated with CL 8 
from t = 0 to t = 5 the CPLI time is shown in Table 6 . The results are roughly consistent with an 
eighth-power dependence of the accuracy on the time step At, expected for CL 8 . 

Let us finally discuss the durations of the runs for the extreme resolution of 8192 2 harmonics, 
used to calculate the radii 5 of the cylinder of analyticity, shown in Fig. 5. Since in a Lagrangian 
code the time step does not depend on the spatial resolution, we expect even larger CPU gains than 
reported above. This is indeed the case: at the resolution of 8192 2 harmonics and in double precision, 
the integration from t — 0 to t — 5 with the four-mode initial condition and our CL 8 algorithm took 
5.5 hours of CPU on the computer mentioned at the beginning of this section, whereas on the same 
computer, the standard RK4 runs for 564 hours, a little over one hundred time slower. The vorticity 
discrepancies between the two outputs are less than 10 -13 up to t — 4.5, as expected from just 
rounding errors, and increase to 4.5 x 10 ” 12 at t — 5, which is consistent with the spatial truncation 
errors estimated in Section 4.1. Note that the truncation order S = 8 , used for the Lagrangian 
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Figure 9: Fourier truncation artefacts (tygers) in Eulerian simulations for the 4-mode initial condition with 1024 2 
harmonics by RK4. Isolines of the vorticity (step 0.5) at t = 5 (a) and of its Laplacian (step 500) at t = 4 (b). 
Blow-ups of these figures: vorticity (step 0.05) at t = 5 (c) and its Laplacian (step 0.5) at t = 4 (d). 

time-Taylor series is well below the optimal order, which is around S' = 20. Using another computer 
with a larger memory, allowing to run CL20, we timed CL20 and RK4 over identical time intervals 
and found that the former runs about 200 faster than the latter. Of course the exact ratios depend 
not just on memory size, but also on the presence or not of high-wavenumber modes in the initial 
condition and on the precise method of parallelisation. Such issues will be discussed elsewhere. 

4 . 4 . Galerkin truncation artefacts 

In spectral simulations of ideal flow, one typically runs out of spatial resolution after some time 
of integration. As pointed out in Section 3.2, for analytic initial data, spatial truncation errors 
roughly behave as , which ceases to be small when complex-space singularities approach the 

real domain to within a few spatial grid mesh sizes. In Eulerian simulations, localised artefacts, in 
the form of spurious oscillations of contour lines, called “tygers,” are then appearing in unexpected 
locations [39]. These are produced by truncation-generated waves interacting resonantly at places 
where the fluid velocity matches that of developing fine-scale structures. For the 4-mode initial 
condition, the phenomenon is illustrated in Fig. 9, which shows Eulerian RK4 simulations with the 
resolution of 1024 2 harmonics at the earliest times when tygers become visible. As is typical for 
small-scale objects, the tygers become at first conspicuous in contour levels of the Laplacian of the 
vorticity, around t = 4 when 6 is about 5 spatial mesh sizes, before they are seen in the vorticity 
itself, around t = 5 when 5 is about 1.5 spatial mesh sizes. For the random initial condition at the 
same resolution (not shown) similar phenomena happen around t = 0.71 and t — 1. 

Do analogous artefacts emerge in our Cauchy-Lagrangian simulations? Contours of the Laplacian 
of the vorticity show indeed similar spurious oscillations (for the same flow and the same spatial 
resolution) around t — 4.1 (see the last panel in Fig. 8). However, the vorticity shows nothing 
comparable until 5 drops well below one spatial mesh. 

Hence, Galerkin truncation artefacts are sharply reduced in our CL simulations, compared to 
Eulerian ones with the same resolution. This is surprising given that the Eulerian and Lagrangian 
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Figure 10: L 2 norms of coefficients of the displacement, ||£ < ' s - ) ||, and vorticity, ||w s ||, in the Lagrangian (a) and Eu- 
lerian (b) time-Taylor series, respectively, computed at t = 0 for the 4-mode initial condition. Solid line: double 
precision, dashed line: quadruple precision. Thin lines: resolution 512 2 harmonics, thick lines: 1024 2 harmonics. 

coordinates are resynchronised at each time step, so that not much difference in behaviour and 
deterioration in time of solutions would be expected. Part of this reduction could be due to the 
smoothing effect of the interpolation stage, also observed in the semi-Lagrangian approach of [32], 
There is, however, another possibility: solutions truncated to a finite number of spatial Fourier 
modes (generally referred to as Galerkin truncated) behave differently in Eulerian and Lagrangian 
coordinates. In Eulerian coordinates, it is easily shown that quadratic invariants (such as energy and 
enstrophy, in 2D, or energy and hclicity, in 3D) remain invariant after Galerkin truncation. Hence, 
as observed in [39], energy and enstrophy, which in the absence of truncation would be transferred to 
smaller scales, have to be transferred somewhere else and get pushed into the tygers. In Lagrangian 
coordinates, the behaviour is different: The energy is E — (1/2) f \x(a, t)\ 2 d 2 a, where the integral 
is over the [0, 27 t] x [0, 27t] periodicity domain. Before truncation, the energy is exactly conserved in 
any time interval At. After Galerkin truncation, it is only conserved up to second order in At, and 
the same holds for enstrophy. 

5. Rounding noise in high-order Lagrangian vs Eulerian time-marching 

Some challenging open problems in hydrodynamics, such as the issue of blow-up for 3D flow [22] 
or for axisymmetric 2D flow [31], may require extremely accurate simulations in order to distinguish 
genuine blow-up from simulation artefacts. This means that spatial and temporal truncation errors 
must be kept very small and the precision very high. 

For flow in a domain of sufficiently simple geometry (foremost space-periodic flow), the most ac¬ 
curate methods for handling the spatial discretisation are of the spectral or pseudospectral type [23]. 
The reason is that, for spatially analytic flow, the Fourier coefficients fall off exponentially at large 
wave numbers [17, 41] and, thus, the truncation error drops exponentially when increasing the res¬ 
olution. Unfortunately, (pseudo-)spectral simulations hardly ever use temporal discretisation (time- 
marching) of order higher than 4, the most frequent being Runge-Kutta of fourth order (RK4). We 
have therefore a striking imbalance between temporal and spatial discretisation. Of course, when 
using Eulerian methods, one generally has to satisfy the Courant -Friedrichs-Lewy [12] condition: 

k m ax U rn ax A t < C s , (26) 

where L max is the maximum wave number (with dealiasing taken into account), I/ max is the maximum 
velocity in physical space, and C s is an order unity constant, depending on the numerical scheme. As 
a consequence, high spatial resolution requires a very small time step, so that a fourth-order scheme 
might be deemed sufficiently accurate. 

The situation is drastically different when working in Lagrangian coordinates. The time step is 
not controlled by a CFL condition any more, but solely by the radius of convergence of the timc- 
Taylor expansion for the Lagrangian map. This radius is roughly the inverse of the largest velocity 
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gradient present and does not depend on spatial resolution. Hence, in principle, there is no reason 
to limit oneself to low-order time-marching. 

We now report results of experimenting with very high-order time-Taylor expansions in both 
Lagrangian and Eulerian coordinates. We will see that the severest limitation stems from the growth 
of rounding errors. 

In Fig. 10, the L 2 norms of time-Taylor coefficients in Lagrangian coordinates, obtained by the 
Cauchy-Lagrangian method, and in Eulerian coordinates, obtained by the Eulerian Taylor method, 
are shown for the 4-mode initial condition, computed in double and quadruple precision up to the 
orders beyond which an unnatural behaviour, dominated by rounding errors, sets in. Taylor coef¬ 
ficients of the displacement, ^ s \ and of the vorticity, u> s , in the expansions of the Lagrangian and 
Eulerian solutions, respectively, have been computed around t = 0 with two resolutions of 512 2 and 
1024 2 harmonics. 

The most striking feature is that the L 2 norms show a range of orders over which the results 
display very little dependence on precision and then a rapid transition to a regime in which drastic 
differences of many orders of magnitude are observed. This happens around a “transition order” 
that depends very much on precision and which is considerably larger (by about a factor 4) in the 
Lagrangian case than in the Eulerian case. (There is also some dependence on the resolution.) The 
exponential behaviour (straight line in the lin-log coordinate system used) is not observed at the 
smallest values of the order s. Indeed, such exponential behaviour ~ R~ s of the coefficient of order 
s of a power series with radius of convergence R is typically only an asymptotic property at large s. 

The large difference between the Lagrangian and the Eulerian computations in the transition 
order, where rounding noise takes over, can be seen as a consequence of a phenomenon discovered 
by Ebin and Marsden [14], which is called the loss of derivatives in the Eulerian formulation. In this 
formulation, each time derivative dt is accompanied by a space derivative v ■ V. As a consequence, if 
the initial velocity has only a finite number of space derivatives, it will also have a finite (generally, the 
same) number of time derivatives and, thus, cannot be analytic in time. In Lagrangian coordinates, 
this is not the case: time-analyticity of the Lagrangian map, for at least a finite time, holds if the 
initial velocity is slightly smoother than once differentiable (see Section 2). Ebin and Marsden, 
although they did not prove time-analyticity in the Lagrangian formulation, did show that the 
variational formulation of the Euler equation, which is essentially Lagrangian, allows interpreting it 
as an ordinary differential equation in function spaces, such as Sobolev spaces, possessing a finite 
number of space derivatives. The proof was given using the language of differential geometry on 
infinite-dimensional Riemannian manifolds. Although somewhat simplified by Bourguignon and 
Brezis [4], the proof did not suggest any numerical implementation. This became possible with the 
rediscovery of Cauchy’s Lagrangian formulation, enabling us to generate time-Taylor coefficients of 
arbitrary orders. 

Now, we comment on why loss of derivatives in Eulerian coordinates is responsible for rounding 
noise that grows very fast with the order of the time-Taylor coefficient. The 4-mode initial condition 
used for the present study has of course space derivatives of arbitrary order (it is an entire function). 
Hence time derivatives of arbitrary orders at t — 0 also exist. It is easy to see that spatial Fourier 
harmonics of the time-Taylor coefficient of order s (not subjected to Galerkin truncation) vanish 
beyond a wave number growing linearly with s. (This is true both in Lagrangian and Eulerian 
coordinates.) Once the problem is discretised in space and FFTs are used, those wave vectors which 
should not be excited will be populated by rounding noise. In the Eulerian case, this noise is amplified 
by the presence of the space derivatives in the recurrence relations, each yielding a factor i k on the 
Fourier coefficients. In the Lagrangian case, the recurrence relation is written in terms of gradients of 
time-Taylor coefficients, but otherwise we only perform additions, multiplications and apply bounded 
Calderon-Zygmund operators. Thus, amplification of noise at high wave numbers is not so strong 
as in the Eulerian recurrence relations. A further confirmation of this mechanism is obtained by 
considering the relative change with the resolution of the order s at which the transition to noise 
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Figure 11: Same quantities as in Figs. 10 (a) and (b), but computed for time-Taylor expansions around various times 
(as labelled). 


takes place. When changing from resolution of 512 2 to 1024 2 harmonics, the transition orders change 
by just a few percent in the Lagrangian calculations and by nearly 20% in the Eulerian calculations. 
(Of course, the higher the resolution, the more operations and thus the more accumulation of rounding 
errors takes place.) 

The advantage of Lagrangian time-marching over Eulerian time-marching from the point of view 
of rounding noise is not limited to initial conditions comprising a finite number of non-vanishing 
harmonics. To show this, we use a multistep method with the same initial data as above, but where 
a very broad band of Fourier harmonics become excited beyond the first time step. We begin by 
explaining how we choose our time steps in the Eulerian calculations. The initial data considered 
here are analytic in the space variables and remain analytic in space and time as long as there is 
no blow-up (for ever in 2D) [3]. At any time t > 0 there is a finite radius of convergence R(t) of 
the Eulerian time-Taylor series giving the solution (say, the vorticity) at time t + At in terms of 
the solution at time t. Certainly, we have to take At < R(t), but this is not enough, because of a 
rounding noise problem, closely linked to what we discussed above. 

In Eulerian coordinates, the fastest temporal variation of small-scale eddies comes from them be¬ 
ing swept by the large-scale energy-carrying eddies. Denoting by tt{t) the amplitude (say, the vortic¬ 
ity) associated with such an eddy, the sweeping can be modelled by the equation d t d + ( U ■ V)d = 0, 
where U is a large-scale velocity which can be taken uniform to leading order (the expansion pa¬ 
rameter being the ratio of scales of the swept eddies and the sweeping eddies). Let us denote 
by dfc the Fourier coefficients of the small-scale vorticity. The sweeping equation has the solution 
dfc(f + At) = exp(—ifc ■ U At)ttk(t). Expanding the exponential in a Taylor series in At, we find: 

4(t + Ai) = 4(t)f: Hfc '( fAt) ‘ . (27) 

s=0 S - 

We define C(k , U) = | k ■ U Af|. Its maximum over all wave vectors and velocities present in the 
flow is the Courant number Co = k max U ma , x At. Here, /c max is the maximum wave number, determined 
by the number of grid points per spatial period, N, and the dealiasing (typically, /c max = N/ 3), and 
Cmax is the maximum velocity. We observe that, in the series (27), the moduli of the various terms 
decrease with s when Co < 1. Otherwise, for sufficiently small eddies, swept by sufficiently fast large 
eddies, i.e., for C{k,U) > 1, the moduli increase at first with the order s while s < C(k,U). Of 
course, the whole sum, being an imaginary exponential, has a unit modulus, but this is a sum of 
a large number of terms, some of which have moduli that are exponentially large in the Courant 
number. The accuracy of such a summation for large Co is low, due to the finite precision of the 
computations. This can result in catastrophic loss of precision - a well-known phenomenon that can 
mar numerical summation of real series with alternating signs (see, e.g., [1]). To avoid this difficulty, 
one can require Co < 1, which is precisely the standard Courant-Friedrichs-Lewy condition [12], 
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Note that the CFL criterion for the Eulerian Taylor method does not stem from a stability condition, 
but from the finite precision of the computations. 

Figure 11 shows the L 2 norms of time-Taylor coefficients of various orders, in Lagrangian and Eulc- 
rian coordinates, respectively. All computations are done in quadruple precision with the resolution of 
1024 2 Fourier harmonics. The various graphs are plotted for Taylor expansions around different out¬ 
put times. Because of the constraints on time steps that are quite different in Lagrangian and Eulerian 
coordinates, the number of time steps between successive output times is quite different in Lagrangian 
coordinates (typically of the order of tens) and in Eulerian coordinates (typically of the order of hun¬ 
dreds). While the Lagrangian and Eulerian time-stepping is performed using Taylor expansions 
truncated at the eighth order, at output times they are pushed to much higher orders. For the Eu¬ 
lerian integrations we observe again a transition to rounding-noise dominated coefficients at orders 
roughly in the range 15-20 (a little beyond order 20 for t — 3). In the Lagrangian case this also 
happens, but at orders about four times larger (not shown in the figure). Thus, we observe essentially 
the same phenomenon as for the Taylor expansions around t — 0. The explanation is now somewhat 
different: in the Eulerian case it is not just the noise that is amplified by the i k factors stemming 
from spatial differentiation, but the actual signal, which can now be much stronger than the direct 
rounding noise. However, when we use FFTs, there is indirect rounding noise, stemming from much 
lower wave numbers that have much higher amplitude. Furthermore, in the Eulerian simulations, this 
indirect rounding noise gets again amplified by the i k factors stemming from the space derivatives. 

6. The radius of convergence and depletion 

A well-known advantage of Lagrangian and semi-Lagrangian methods is that the time step is 
constrained not by the usual CFL condition but, roughly, by the inverse of the largest velocity 
gradient. For our Cauchy-Lagrangian method, the optimal time step At for a given accuracy, when 
starting from time t, is a fixed fraction of the radius of convergence Ri n f(t) of the time-Taylor series 
for the Lagrangian map (see Section 3.5). It is of interest to know how A? in f(t) changes in time, since 
this determines how the time step should be changed and also the total number of time steps that 
will be needed to reach a given termination time tfi n . When using initial conditions that have mostly 
low-wave-number harmonics, such as the 4-mode flow, the subsequent transfer of excitation to higher 
harmonics will be accompanied by some growth in the velocity gradient: substantially so in 3D and 
more moderately in 2D, because of the conservation of vorticity. 

As a consequence, we expect a decrease of the radius i? in f(t). If i? in f(t) —> 0 as t —> from 
below, a blow-up of the solution takes place at the finite time £*. Whether the converse holds is so 
far unknown. Anyway, in the present numerical studies, we are only concerned with the 2D case. 
Before turning to numerical results, let us recall the proven constraints. Wolibner’s theorem [42] 
then guarantees regularity for all t > 0 when the initial vorticity is Holder-continuous of exponent 
a 0 > 0. More precisely, Wolibner showed that, at time t > 0, the vorticity is Holder-continuous of 
exponent a(t ) > oiq e~* sup ^°l, where sup |cao| is the supremum of the modulus of the initial vorticity. 
The constraint on vorticity conservation does not prevent the velocity gradient from growing, but it 
cannot grow faster than e * sup , a bound which is actually sharp [27]. Consequently, for such 2D 
flow, the radius of convergence has a lower bound proportional to e —* sup l^ol ( see Section 2). 

The actual numerical results, shown in Fig. 12, indicate that the decrease in time of the radius of 
convergence is much slower than this bound. Consider, for instance, the 4-mode initial condition, 
for which the supremum of the initial vorticity is 2.8. The aforementioned lower bound could have 
the radius of convergence decreasing by a factor e~ 5x2 ' 8 « 8 x ICC' over the time span from 0 to 5. 
Actually, the radius of convergence decreases only by about 30%. As to the random initial condition, 
its radius does not decrease at all. For comparison, the radii of convergence for the Eulerian Taylor 
method are also shown: they decrease much faster. 

This huge discrepancy between the proven bounds and the actual results is most likely related to 
the phenomenon of depletion: in both 2D and 3D incompressible flow it is frequently found that the 
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Figure 12: Temporal evolution of the radii of convergence Ri n {(t) of the time-Taylor series of the displacement £(a,t) 
simulated by the Cauchy-Lagrangian method (dashed lines), and of the vorticity u>(a,t) simulated by the Eulerian 
Taylor method (solid lines). Dots show the times at which the radii were computed. Initial condition: 4-mode (a), 
random (b). Resolution: 1024 2 harmonics, quadruple precision computations. 


nonlinear effects are growing in time much slower than permitted by the best proven bounds [19], 
probably because such flow tends to organise itself into structures that have an almost vanishing 
nonlinearity, e.g., quasi-one-dimensional flow depending mostly on a single coordinate. Depletion is 
definitely why 2D flow is often found numerically to be even more regular than implied by Wolibner’s 
lower bound, and also why 3D simulated flows show little reliable evidence for blow-up (one possible 
exception is a recently studied swirling axisymmetric flow with a solid boundary [31]). 

7. Conclusion 

We have developed a novel Cauchy-Lagrangian (CL) method to solve the initial-value problem 
for ideal incompressible flow, governed by the Euler equation. The method operates in Lagrangian 
coordinates, relies on temporal analyticity of particle trajectories in time and employs the recurrence 
relation (9) for coefficients of the Lagrangian time-Taylor series, stemming from the Cauchy invari¬ 
ants. Our experiments in 2D show that this method is highly efficient (like many semi-Lagrangian 
methods, it does not suffer from the Courant-Friedrichs-Lewy constraint) and enables typically an 
order of magnitude faster time-marching than the fourth-order Runge-Kutta method and the Eulc- 
rian Taylor method, at the same time being significantly more accurate. As shown in Section 4.3, 
the gain in speed can even be higher, in the range 100-200, at the highest spatial resolutions. Last 
but not least, the Cauchy-Lagrangian method allows very-high-order time marching without the 
severe rounding error problems that affect high-order Eulerian time marching. Of course, these im¬ 
pressive results are so far obtained for just a narrow niche of problems: at present, incompressible 
space-periodic two-dimensional flow. 

The CL algorithm has been presented in Section 3 for flow of arbitrary dimension. Since the 
CL method allows efficient high-order and high-precision simulations with only modest growth of 
rounding errors, it is a strong candidate for the reliable numerical exploration of the blow-up problem 
for the 3D Euler equation (see, e.g., [22]). 

The method, here described for incompressible Euler flow, can be readily extended for solving 
certain other problems arising in various areas of physics. These are described by equations for which 
known theorems guarantee the time-analyticity of Lagrangian fluid particle trajectories. Analytic¬ 
ity in time of the Lagrangian map for at least a finite time, under the assumption that the initial 
velocity gradient is, e.g., Holder-continuous, holds for the Euler equation (which applies to inviscid 
incompressible fluid flow and to the 2D drift-Poisson equation for an electron plasma). Furthermore, 
it also holds for the surface quasi-geostrophic equation, the incompressible porous medium equa¬ 
tion and the Boussinesq equation. So far, the latter three have been handled not by a recurrence 
relation approach, using Cauchy’s Lagrangian formulation, but by applying the Faa di Bruno for- 
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mula [11] which gives an explicit but somewhat cumbersome representation of all the time-Taylor 
coefficients [10]. Analyticity results, derived from recurrence relations similar to (7)-(8), were also 
proved for compressible flow of cosmological interest, such as the Euler-Poisson equations for an 
Einstein-de Sitter universe [44], recently extended to the more realistic case of a ACDM universe 
[38]. The Euler-Poisson equations apply to the collisionlcss flow of dark matter only as long as 
multi-streaming (development of caustics with more than one velocity) is absent or negligible. In 
such a framework, our methods are useful in at least two cases: the blow-up problem (called shell¬ 
crossing by cosmologists) and the reconstruction problem, in which the whole dynamical history of 
the universe is obtained from just the spatial distribution of matter at the present epoch [18, 7]. 

But a warning is necessary: not all ideal fluid flow problems are endowed with time-analytic 
Lagrangian fluid particle trajectories. For example, this is apparently not the case of magnetohydro¬ 
dynamic (MHD) ideal flow and of barotropic flow. For non-diffusive MHD flow, one can still write 
recurrence relations for the time-Taylor coefficients of the Lagrangian map. However, the Lagrangian 
gradient of a time-Taylor coefficient of a higher index does not involve just the gradients of those of 
lower index, but also their space derivatives of order up to three. In other words, in contrast with the 
ordinary Euler equation which does not lose derivatives in Lagrangian coordinates (see Section 5), 
the MHD equations do, and thus MHD flow with only limited initial spatial smoothness lacks ana¬ 
lyticity of fluid particle trajectories. Such a loss of derivatives may be viewed physically as due to the 
frozenness of the magnetic field. The barotropic equations govern compressible fluids, in which the 
pressure is a function of density only. They play an important role in geophysical fluid dynamics [33]. 
The standard derivation of the Cauchy invariants applies also to the barotropic equations. However, 
the second equation, stemming from mass conservation, is profoundly modified and, in it, loss of 
derivatives and hence of analyticity occurs again. 

Another class of flows which might be amenable to some variant of the Cauchy-Lagrangian 
method is viscous flow (e.g., Navier-Stokes flow) that is initially analytic in the space variable (if it 
is not space-analytic initially, it becomes so at an arbitrary short positive time [16]). Such flow is 
also analytic in time in both Eulerian and Lagrangian coordinates. There is however a new practical 
problem: when the equations are written in Lagrangian coordinates, if we try following Cauchy’s 
approach, we find that viscous terms involve nonlinearities in the Lagrangian map of degree four in 2D 
and of degree 6 in 3D; thus the recurrence relations become somewhat cumbersome. An alternative 
and much simpler strategy for viscous flow is a fractional step method [43] in which the reversion 
to Eulerian coordinates by interpolation is followed by a solution of the viscous heat equation. This 
is related to the Trotter-formula approach to Navier-Stokes flow discussed in §13 of [14]. At high 
Reynolds numbers, this method will result in a low-order scheme for the smallest viscosity-dominated 
scales but the larger scales will still benefit from the high accuracy permitted by the CL method. 

In this first paper on the Cauchy-Lagrangian method we have concentrated deliberately on ex¬ 
plaining and validating the technique. We expect that, over the years, this technique will be ap¬ 
plied and improved for a range of situations and problems. It is of course desirable to handle flow 
with boundaries: for very smooth boundaries, this can in principle be done by solving appropri¬ 
ate Hclmholtz-Hodge problems with boundaries. As explained in Section 3, the basic ideas for the 
Cauchy-Lagrangian method are the same for 2D and 3D flow. At high resolution, it becomes impor¬ 
tant to have an efficient parallelisation of both the Lagrangian time-stepping and the interpolation 
phase back to Eulerian coordinates; this is being currently explored. One of the very challenging 
applications of such 3D Euler computations is the possible persistence for all times of an initially 
assumed smoothness vs. the loss thereof after a finite time (blow-up problem). In this context, we 
finally observe that real-space blow-up cannot take place in bounded planar 2D flow, but there is 
already strong evidence that in 2D the Lagrangian time-Taylor series (around t = 0) has a finite 
radius of convergence. Hence, Lagrangian trajectories should have complex-time singularities, whose 
nature can be determined by careful simulations. 
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Appendix A. Bounds for the space-dependent radii of convergence 

The Lagrangian time-Taylor series (6) is at the core of the Cauchy-Lagrangian numerical method. 
For various reasons explained in the body of this paper, it is of interest to know its radius of 
convergence. More precisely, we want to know how to determine the infimum over the Lagrangian 
space of the radii of convergence of the time-Taylor series. A theorem, announced in Section 2.2, 
allows us to relate this radius to that of the ordinary Taylor series obtained by replacing each term 
in (6) by its L p norm. We give now a formal statement and proof of this theorem. 

We consider an abstract infinite function series 

OO 

Ua ,t) = Y,C ( ‘\a)t‘, (A.l) 

S = 1 

where a is in the d-dimensional periodicity domain (torus) T d . This need not be a Lagrangian 
solution of the Euler equation; in particular, it suffices for our purposes to provide arguments for a 
scalar-valued series. Let R(a) denote the radius of convergence of the time series at point a. We 
define i? in f = essinf a i?(a); we denote the L p norm of a function £(a) by ||£(a)|| p and the radius of 
convergence of the series ||£^(a)||p£ s by R p . The Lebesgue space L p (T d ) under consideration 
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is that of scalar-valued functions of the spatial variable a of an arbitrary dimension d; we assume 
that all £^( a ) and the sum f(a,t) are space-periodic. 

We prove now the following 
Theorem. For some p > 1, let all £( s \a ) be in If. Then 

1. For all t < R p , the sum £(a,t) belongs to L p , and i? in f > R p . 

2. If £(a,t) is also in If for any t < R, then R p > 

Proof. The radius of convergence of the series (A.l), R(a), does not change, when the coefficients 
£( s )( a ) are replaced by their absolute values, |£^(a)|. Similarly, applying the absolute value does 
not modify the Lebesgue space norms. Consequently, without any loss of generality, we assume 
£( s )(a) > 0 for all s and a, otherwise replacing the coefficients by their absolute values. 

To prove statement 1 of the theorem, we consider the truncated series £s(a,f) — X)f=i ^ s \ a )t s 
for a fixed t in the interval 0 < t < R p . Since the partial sums £5 (a, t) constitute a Cauchy sequence 
in L p , it converges in L p to a limit function £(a,t) also belonging to L p . By a well-known theorem 
of functional analysis (see, e.g., [ 2 ]), one can choose a subsequence fs k (a,t) that converges almost 
everywhere to £(a, t). (In other words, for our fixed t, there exists a set of measure zero such that, at 
any point a in the complement of this set, lim^^ £g fc (a, t ) = £(a,£), the limit being here understood 
as the usual limit for infinite sequences of real numbers.) However, since f( s \a) > 0 for all s > 1, the 
partial sums fs( a ,t) grow monotonically at each point a (for the fixed t) on increasing the index S, 
and hence the entire sequence £s(a, t) converges for S —> oo almost everywhere to £(a,t). Therefore, 
any fixed time t < R p is within the disk of convergence of the series (A.l), i.e., R(a) > R p for almost 
all a, which implies the first statement of the theorem. 

We prove now statement 2. For any p > 1 and any quantities a s > 0 such that a s = 1, 

clearly, a s — 1- Taking a s = (a)t s /f(a, t) (recall ^ s \a) > 0 for all s), we thus find 

^((M(a)tr <(”(“, t) (A.2) 

S>1 

for almost all a and any t such that 0 < t < R(a), and hence for any 0 < t < i?; n f. Inte¬ 
grating this inequality over the periodicity domain yields ^) s>1 ||^ s ^ {a)\\ p p t sp < ||C(a,f)|| p . For 

t > R p = (lim sup^oo ||e (s) (a)||p /s )- 1 , the series in the l.h.s. of the latter inequality diverges, but 
it converges for all t < i? in f, for which the inequality holds. Consequently, i? in f < R p . This concludes 
the proof of the theorem. 

It remains to comment on the significance of the condition in statement 2 that £(a,t) belongs to 
L p for all t < R. One can indeed construct an abstract counterexample showing that the radii are 
not necessarily equal, if this condition is not satisfied. 

Counterexample. For simplicity, we consider functions defined in the interval [0,1]. Suppose 
1 = a 0 > Gq > ... > a s > a s+ i > ... > 0 are real numbers such that linq^oo a s = 0. For a > 1/p and 
s > 1 , consider functions ^ s ^( a ), which vanish outside the interval a s _i > a > a s , and £^(a) = a~ a in 
this interval. Evidently, the series (A.l) converges for any t. i.e., R.- m f = oo. The sum (A.l) belongs to 
the Lebesgue space L p ([0,1]) only for \t\ < R p , since by construction ||^(a, t)\\ p = ^) s>1 ||£^(a)|| p f sp . 
By choosing a suitable sequence of nodes a s , we generate functions £^( a ) of prescribed norms 
ll^^( a )llp = 7 s- Since ||'C ( ' s ' ) (a)||p = (pa — 1 )~ 1 ^ p (a l ~ pa — a^If*) 1 ^, the condition lim^ooas = 0 is 
satisfied whenever X} s >i7s = °°j otherwise, > 0 are arbitrary. For instance, % = 7 s is permissible 
for any 7 > 1. We have thus constructed an example of series (A.l), non-summable in L p ([0,1]), for 
which R,\„f = 00 , while the accompanying series (12) constructed of If norms possesses an arbitrary 
prescribed finite radius of convergence R p = I/7 < 1; for such a series, only the first statement of 
the theorem holds. 

Remark. The Lebesgue space L p (T d ) in the statement of the Theorem can be replaced by a Sobolev 
function space W p (T d ) for p > 1 and q > 0. Its norm is defined as ||£(a)|| 9 ,p = ||(/ — V 2 ) 9 / 2 £( a )|| p; 
where I is the identity operator and V 2 the Laplacian (the spaces for p — 2 are often denoted by 
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H 9 (T d )). The quantity R p is then replaced in the modified theorem by R q , p , the radius of convergence 

of the series ||£^(cOllg,p^ i > and instead of R mi we use R qpn i, the essential infimum over space 

(i.e., over the variable a) and over all q', such that 0 < q' < q, of the radii of convergence of the 
series 

OO 

J2(-\7 2 ) q ' /2 ^ s \a)t s . (A.3) 

S= 1 

Convergence of the series implies convergence of (A.3) in L p for all q' < q and, 

by the arguments in the proof of the Theorem, the inequality R g> ; n f > R q , p - Conversely, convergence 
of (A.l) in W p implies convergence of (A.3) in L p for q' = 0 and q' = q, and by the second part of the 

Theorem the radii of convergence of the two series, HC^HpC, and II ( — V 2 )' 3 ' /2 £ ( ^||pf s , are 

both bounded from below by R qpn f, whereby R qpn f < R q , P - Note that R q>p and R qpn { are monotonically 
decreasing functions of q for any p > 1. The inequality R qtP < R q \ p for some q > q' signals existence 
of a ring within the disk of convergence in W^,, where the series (A.l) does not converge in W p . An 
intriguing problem is to establish whether for a solution of the Euler equation R qpn f = R q > pn f and 
R q , P = Rq', P hold for all sufficiently small q ^ q'. 


Appendix B. Computation of the radius of convergence of a power series 
with positive coefficients 

The implementation of the Cauchy-Lagrangian method does not require the evaluation of the 
radius of convergence of the time-Taylor series employed. However, to test the method, to analyse 
its performance and interpret some of the results (Sections 4-6), we need to compute the radii of 
convergence of the time-Taylor series for both CL and ET algorithms. For estimating these, using the 
result of Appendix A, we can replace the full functional Taylor series by the ordinary 

power series ||£ < A(a)|| t s , where || • || denotes the L 2 norm. To obtain reasonable precision on 

the radii, say, a few percent, this computation actually requires the use of quadruple precision and 
the knowledge of enough Taylor coefficients, usually more than is required for the time-stepping. 
Due to this added complexity, we do not compute the radii at every time step, but just frequently 
enough to be able to monitor how the radii evolve in time, e.g., every tenth time step. 

In this Appendix we describe how we calculate the radii for such power series. 

Here is the problem that we want to solve: 

Given a finite number, S, of coefficients of the power series 

OO 

/W = (B.l) 

S=1 

where f s > 0, we want to estimate its radius of convergence, which we denote by R. For our purposes, 
an accuracy of a few percent suffices. 

To begin with, we recall some classical results. By the Cauchy-Hadamard theorem, 

= lini sup a/|7J. (B.2) 

S—>• OO 


A fast way of computing the radius of convergence is the ratio formula 


1 

R 


lim 

S—>• OO 



(B.3) 


when the limit exists. 

With less than a hundred coefficients of the time-Taylor series at hand, we can obtain the radius 
of convergence with a reasonable accuracy using such formulae only, if we know several terms in the 





Figure B.13: The quantity F s /s (solid line) and its least-square fit (dashed line, see (B.5)) (a). Scaled discrepancy d s 
(B.6) (b). 


large-s asymptotic expansion of the coefficients f s . We should derive the asymptotics of f s = ||£ ( - s )(a)|| 
from the recurrence relation (9) for . Unfortunately, no such results are available so far. To 
proceed, we look at our problem from an entirely different angle and ask ourselves the question: Does 
the analytic function defined by the series j ||£; ( '^(a)|| t s have the “generic” simple structure near 
the disk of convergence, such as many analytic functions (B.l) with real positive coefficients f s have? 

A standard result is that the radius of convergence of a Taylor series is the distance in the complex 
f-plane to the nearest singularity(ies). If there is a single such singularity and the Taylor coefficients 
are all positive, the singularity is on the real positive axis at f* = R. The simplest case is when the 
singularity has a simple structure, e.g., a pole or a power-law branch point that can be characterised 
by a local behaviour oc (t* — t) p , where the exponent p is a negative real number. It may then be 
shown that the Taylor coefficients of large orders s have the following behaviour: 

f s ^'ys a e l3s for s —>■ oo, (B.4) 

where /3 — — In R and a = — p — 1. (The same kind of asymptotics applies to Fourier series; see, e.g., 
[ 8 ].) In such a case, if we are able to extract the leading and subleading behaviour of f s , we may 
then infer not only the radius of convergence, but also the type of the nearest singularity on the real 
axis. One method is the Domb-Sykes plot of f s /f s -1 versus 1/s, where the ^/-intercept gives 1/R 
and the slope near the origin gives the exponent a [13, 24], 

We have applied all the standard methods, recalled above. We found that, when using 50-80 
quadruple-precision Taylor coefficients from our L 2 series, such methods allow a determination of the 
radius of convergence with an error in the range 10%-20%. We also tried more advanced extrapolation 
techniques, such as convergence acceleration methods (see, e.g., [35, 34] and references therein) and 
the asymptotic extrapolation technique [25, 35, 34], which have the potential of capturing both 
leading and subleading asymptotic behaviour. Unfortunately, all these techniques failed to reach the 
asymptotic regimes and thus gave no improvement over the more standard ones. 

We found that for the L 2 series with 50-80 terms a rather straightforward method apparently 
gives a better accuracy. We again assume that the L 2 series has a single singularity and that (B.4) 
holds asymptotically for large s. It follows that the logarithms of the Taylor coefficients, F s = ln/ s , 
are given by 

F s ~ In y T o In s + /3s for s —* oo. (B.5) 

We can determine the coefficients a, [3 and 7 by a least-square fit over the available Taylor coefficients. 
That is, we minimise ^ s=1 d 2 s for the discrepancies 

d s = c + a In s + bs — F s . (B. 6 ) 

We now give an example of such determination for the four-mode initial condition (20). Here 
f s are the L 2 norms of the first 80 time-Taylor coefficients £-b for the displacement, computed at 
t — 0. Figure B.13 (a) shows the values F s /s and the fitted curve (c + a Ins + bs)/s. (Note that 
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F s /s ~ /? = — Ini? for large s.) The least square fit gives a = —1.94, b = —0.187 and c = 0.116. 
Thus, the radius of convergence is R = e~ b ~ 1.2 . To estimate the error on the radius R, we show in 
Fig. B.13 (b) the scaled discrepancy d s /s for large s. We see that, beyond s = 20, the fluctuations 
do not exceed 1%. Similar determinations of the radius R(t) for the same initial four-mode flow and 
0 < t < 5 have been made, using between 80 and 50 Taylor coefficients, yielding errors between 1% 
and 3%. 

We have used the same fitting technique to measure the radius 6 of the cylinder of spatial analyt- 
icity around the real domain (in both Lagrangian and Eulerian coordinates). The analogue of (B.4) 
is then the following expression for the shell-averaged energy spectrum E(K), 

E(K) « CK n e~ 25K . (B.7) 

How can we achieve a much better accuracy on radii of convergence? We need to know a lot 
more than 80 coefficients, say, a few hundred. Unfortunately, as observed in Section 5, in quadruple 
precision, rounding noise becomes explosively large just above order 80. We can go a little further by 
avoiding FFTs and performing the Fourier-space convolutions explicitly, but this is computationally 
very expensive. Otherwise, we have to resort to higher than quadruple precision. Recent algorithmic 
developments on FFTs in “medium” precision (up to 400 bits) offer interesting possibilities [26]. 
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